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In an earlier study, we reported on the excitation of large-scale vortices in Cartesian hydrodynamical convection models 
subject to rapid enough rotation. In that study, the conditions of the onset of the instability were investigated in terms of 
the Reynolds (Re) and Coriolis (Co) numbers in models located at the stellar North pole. In this study, we extend our 
investigation to varying domain sizes, increasing stratification and place the box at different latitudes. The effect of the 
increasing box size is to increase the sizes of the generated structures, so that the principal vortex always fills roughly 
half of the computational domain. The instability becomes stronger in the sense that the temperature anomaly and change 
in the radial velocity are observed to be enhanced. The model with the smallest box size is found to be stable against 
the instability, suggesting that a sufficient scale separation between the convective eddies and the scale of the domain is 
required for the instability to work. The instability can be seen upto the co-latitude of 30 degrees, above which value the 
flow becomes dominated by other types of mean flows. The instability can also be seen in a model with larger stratification. 
Unlike the weakly stratified cases, the temperature anomaly caused by the vortex structures is seen to depend on depth. 
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1 Introduction 

Hydrodynamical Cartesian convection simulations subject 
to high enough rotational influence (Co ^ 3) and exhibit- 
ing large enough Reynolds number (Re 30) have been 
reported to generate vortices, the sizes of which are large 
compared to the size of the convection cells (e.g. Chan 
2003, 2007; Kapyla, Mantere & Hackman, 2011, hereafter 
KMH1 1). In the moderate Coriolis number regime, the vor- 
tices are cyclonic, suppressing the energy transport by con- 
vection, and thereby appearing as cooler than their sur- 
roundings. When Coriolis number is increased even further, 
anticyclonic vortices are preferred, enhancing the convec- 
tive energy transport, making the vortices appear as regions 
warmer than their surroundings (KMH1 1). 

In our previous study (KMH11), we proposed that 
such vortical structures could be responsible for cool/hot 
starspots in rapidly rotating late-type stars possessing outer 
convection zones. We were prompted to look into such a 
possibility due to the decorrelation of the surface tempera- 
ture maps, obtained by Doppler-imaging techniques, from 
the distribution of surface magnetic fields, derived through 
Zeeman-Doppler imaging methods (e.g. Donati & Collier 
Cameron 1997; Donati 1999; Hussain et al. 2000; Jeffers et 
al. 201 1; Kochukhov et al. 201 1). 

The resulting temperature anomaly due to the vortices 
was shown to be of the order of 5 percent, being somewhat 



weaker than the temperature contrasts deduced from obser- 
vations. The model, however, was very simple: for example, 
the density stratification in the radial direction was only of 
the order of 23. Due to the low growth rate of the instability, 
requiring several thousand turnover times to saturate, only 
a very limited parameter range was studied. In this study, 
we extend the previous one by investigating a model with a 
larger stratification, models with varying domain size, and 
place the computational domains at different latitudes. 

2 Model 

Our model is based on that used by Kapyla et al. (2009) and 
KMH 1 1 . A rectangular portion of a star is modeled by a box 
situated at colatitude 6. The box is divided into three layers: 
an upper cooling layer, a convectively unstable layer, and a 
stable overshoot layer (see below). We solve the following 
set of equations for compressible hydrodynamics: 
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where B/Bt = d/dt + U ■ V is the advective time deriva- 
tive, v is the kinematic viscosity, K is the heat conductivity, 
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p is the density, U is the velocity, g = —gz is the gravita- 
tional acceleration, and ft = fio ( — sin 8,0, cos 8) is the ro- 
tation vector. The fluid obeys an ideal gas law p = (7 — 1 ) pe, 
where p and e are the pressure and the internal energy, re- 
spectively, and 7 = cp/cy = 5/3 is the ratio of the spe- 
cific heats at constant pressure and volume, respectively. 
The specific internal energy per unit mass is related to the 
temperature via e = cyT. The rate of the strain tensor S is 
given by 

The last term of Eq. (O describes the cooling at the top of 
the domain. Here t(z) is a cooling time which has a profile 
smoothly connecting the upper cooling layer and the con- 
vectively unstable layer below, where t — > 00. 

The positions of the bottom of the box, bottom and 
top of the convectively unstable layer, and the top of 
the box, respectively, are given by (z±, Z2, Z3, z£) = 
(— 0.85, 0, 1, 1.15)d, where d is the depth of the con- 
vectively unstable layer. In the case of larger stratifica- 
tion (Set C), the corresponding vertical positions read 
{z\, Z2, Zs, Zi) — (— 0.4, 0, 1, 1.1) d, resulting in a vertical 
extent somewhat smaller than in Sets A and D. Initially the 
stratification is piecewise polytropic with polytropic indices 
(mi, 7712,7713) = (3, 1, 1), which leads to a convectively un- 
stable layer above a stable layer at the bottom of the domain. 
In a system set up this way, convection transports roughly 
20 per cent of the total flux. Due to the presence of the cool- 
ing term, a stably stratified isothermal layer forms at the top. 
The standard horizontal extent of the box, Ljj = L x = L y , 
is 4c?; the horizontal domain size is varied from half of this 
to double the size in Set A. In Set C with larger stratifica- 
tion, the horizontal extent of the box is 5d, The simulations 
in Sets A and C are made at the North pole, corresponding 
to 6 = 0°, while in Set D, the co-latitude is varied with 
coarse steps to cover the latitude range down to 8 = 60°. 
The simulations were performed with the PENCIL CODlfl 
which is a high-order finite difference method for solving 
the compressible equations of magnetohydrodynamics. 

2.1 Units and non-dimensional parameters 

Non-dimensional quantities are obtained by setting 

d = g = Pa = c P = 1 , (5) 

where po is the initial density at z% . The units of length, 
time, velocity, density, and entropy are 

[x]=d, [t\ = y/dfi, [U] = ^d~g, 

[p] = Po , M = c P . (6) 

We define the Prandtl number and the Rayleigh number as 

v gd i ( 1 ds\ 

Pr = — , Ra=^— -—T-) , (7) 

where \o = K/(p m cp) is the thermal diffusivity, and p m 
is the density in the middle of the unstable layer, z m — 



i(z3 — z%). The entropy gradient, measured at z„ 
non-convective hydrostatic state, is given by 



in the 
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where V — V a d is the superadiabatic temperature gradient 
withV a(1 = I-I/7, V = (d\nT/d\np) Zin , and where Hp 
is the pressure scale height. The amount of stratification is 
determined by the parameter £0 = (7 — l)eo/(gd), which 
is the pressure scale height at the top of the domain normal- 
ized by the depth of the unstable layer. We use £n = 1/3 in 
Sets A and B, which results in a density contrast of about 
23 across the whole domain, and roughly 9 over the con- 
vectively unstable layer. We make one run with higher strat- 
ification (CI), for which £0 = 1/6, resulting in a density 
contrast of roughly 230 over the convectively unstable layer. 
We define the Reynolds and Peclet numbers via 

^rms t-. ^rms 



Re = 



Pe = 



= Pr Rc 



(9) 



where k{ — 2n/d is adopted as an estimate for the 
wavenumber of the energy-carrying eddies, and u rms = 
y/3u*. This definition of u lms neglects the contributions 
from the large-scale vortices that are generated in the rapid 
rotation regime. Note that with our definitions Re and Pe are 
smaller than the usual ones by a factor of 27r. The amount 
of rotation is quantified by the Coriolis number, defined as 

Co = °- . (10) 

We also quote the value of the Taylor number, 

Ta = (2fl d 2 /is) 2 , (11) 

1/9 

which is related to the Ekman number via Ek = Ta ' . 
2.2 Boundary conditions 

The horizontal boundaries are periodic for all variables. 
Stress-free conditions are used for the velocity at the ver- 
tical boundaries. 

U x , z - U y , z =U Z = 0. (12) 

The temperature is kept constant on the upper boundary and 
the temperature gradient 

dT = -g 
dz cy(7 — l)(m + 1) 



(13) 



http://code.google.eom/p/pencil-code/ 



is held constant at the lower boundary, yielding a constant 
heat flux Fq = —KdT/ dz through the lower boundary. 

3 Results 

In an earlier study, we investigated the excitation of large- 
scale vortices in Cartesian domains with weak density strat- 
ification, located at the North pole of a rapidly rotating 
star (KMH11). In that study, we investigated the limit- 
ing Reynolds and Coriolis numbers (i.e. inverse Rossby 
numbers) above which the instability was excited. For the 
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Table 1 Summary of the runs. Stars indicate that the simulation has not been run to a saturated state. The amount of 
stratification, Ap, is measured over the convectively unstable region. The temperature anomaly, T cyc /T, where T cyc is 
the extremum of temperature in a cyclonic region and T the mean temperature of certain horizontal layer, is measured at 
the middle of the convective layer, z m . The dimensionless input heat flux at the lower boundary of the box is given by 
Fo = Fo/(p c s), where c s is the adiabatic sound speed and p is the density, both measured at the lower boundary of the 
domain. The last column indicates the presence of cyclonic (C), anti-cyclonic (A), or both types (A+C) of vortices. 



present study, we have extended our analysis to models with 
varying computational domain size (Sect. 13. U . place the 
box at different latitudes of the star (Sect. 13.2b . and present 
a model with higher density stratification (Sect. 13.3b . The 
most relevant input parameters for the runs and some essen- 
tial output quantities are listed in TableQ] As is evident from 
this table, some of the runs generating vortices have not yet 
reached a completely saturated state (these runs are marked 
with stars), even though they have been integrated for sev- 
eral thousands of turnover times. In such cases, quantities 
such as the temperature anomaly between the vortex and its 
surroundings, might still be underestimated, and the listed 
numbers in Table Q] should be regarded as lower limits. 

3.1 Dependence on the domain size 

In our previous computations with the standard box size of 
Ln=4:d (KMH1 1), we observed a clear tendency of the sizes 
of the vortices to approach the wavenumber k/k\=l, i.e. 
they tended to fill in as large a fraction of the horizontal ex- 
tent as possible. This prompted us to study the dependence 
of the instability on the horizontal extent of our Cartesian 
box. In Run A3 of this paper (see TableQ]), both the horizon- 
tal extent and resolution of the domain are halved. Run A2 
is a model with the standard box sixe (actually otherwise 
identical to D5 except for the value of the kinematic viscos- 
ity), and in Run A3 the resolution and horizontal extent are 
doubled. In all the runs, the computational domain is located 
at the North pole of the star. 

In Run A3, the Reynolds (roughly 39) and Coriolis num- 
bers (roughly 8.7) are clearly above the critical values found 
in KMH11; still, no vortices are excited. This is evident 
from the leftmost panels of Fig. [TJ where we show the tem- 
perature field (upper panel) and radial velocity (lower panel) 
at the middle of the convectively unstable layer, z m . Some 
large-scale kjk\ = 1 pattern can be detected in the tem- 
perature field, that might be indicative of the early stages 
of the vortex-instability. This run, however, was continued 



up to 6500 turnover times, matching the timepoint of the 
slice plotted in Fig. Q] It is not completely ruled out that a 
very slowly growing vortex instability mode is present, but 
in any case its growth rate is strongly reduced from the stan- 
dard box runs presented in KMH1 1 . 

Run A2 shows very similar behavior to the earlier cal- 
culations presented in KMH11: the Reynolds number is 
clearly supercritical to the instability (roughly 44), and 
Coriolis number (roughly 7.7) in the regime where the ex- 
citation of an anti-cyclone, i.e. a vortex rotating in opposite 
direction to the overall rotation of the domain, was reported. 
As shown in the middle panels of Fig. Q] a vortex rotating in 
the clock-wise direction is seen in the velocity field (lower 
panel), the structure being warmer than its surroundings 
(upper panel). The temperature anomaly is slightly less than 
five percent, very close to the number reported in KMH1 1. 
The temperature across the vortex, normalised to the mean 
temperature of the horisontal layer, is plotted in the upper 
panel of Fig. [2] with a red, dashed line. In the region of the 
anti-cyclone, the vertical velocities become somewhat en- 
hanced, as evident from the lower panel of the same figure, 
where the velocity profile normalised to the rms velocity of 
the horisontal layer across the vortex is plotted with a red, 
dashed line. 

In the model with the largest horizontal extent, Run Al, 
exhibiting very similar Reynolds and Coriolis numbers in 
comparison to the standard box Run A2, two instead of 
one anticyclones are seen. The growth rate of the instabil- 
ity, measured from the growth of the horizontal velocity 
components, increases by a factor of 2.3 from Run A2 to 
Al. As evident from the rightmost panels of Fig.Q] the two 
structures are of unequal strength, possibly suggesting that 
the system has not yet completely saturated, although it has 
been followed for up to nearly 3000 turnover times. Based 
on earlier experience with the standard box runs (KMH1 1), 
where multiple vortices were also seen, when followed long 
enough, only one vortex very close to the k/k\=\ wavenum- 
ber would be the stable end point for this value of Co. 
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Fig. 1 Temperature (upper row) and radial velocity U r (lower row) in the middle of the convectively unstable layer, i.e. 
at z m , for Runs A3, A2 and Al (from left to right). For Run A3, the slice is taken from the time t — 6400r to , where 
T~to = (u lms k{) . For Runs A2 and Al, the snapshots are taken at 3400r to and 2700Tt o . 



The temperature and vertical velocity anomaly across the 
stronger vortex is plotted in Fig. [2] with solid, black lines. 
As can be seen from this figure, the heating effect in the 
middle of the vortex is roughly twice as large as in the stan- 
dard box case, suggesting that the strength of the instability 
is indeed dependent on the box size. No such dramatic dif- 
ference can be seen in the vertical velocity cut across the 
vortex (Fig. |2]lower panel). 

The size of the stronger vortex in Run Al, again, is 
nearly half of the domain size, i.e. approaching the k/k\=l 
mode. In Fig. [3] we plot the power spectra for the kinetic 
energy from Runs Al (black lines) and A2 (red lines); the 
wavenumber scale for Run Al is re-scaled to match the one 
of Run A2. In the early stages, when no vortices are yet 
excited (dashed linestyles), both the power spectra consis- 
tently peak at intermediate wavenumber of roughly fc/fci=7; 
this number reflects the size of the turbulent eddies due to 
convective motions. Due to the appearance of the vortices, 
the energy contained in large scales grows, and dominates 
the flow in the nearly saturated stage. Most of the power is 
seen near the wavenumber k/k\=\. In Run Al, the insta- 
bility would still have 'space' to advance into even lower 



wavenumbers i.e. larger scales; followed up even further, it 
might still do so. 

3.2 Latitudinal dependence 

In Set D, we place the computational domain at different co- 
latitudes with a very coarse latitude grid of [0, 15, 30, 45, 60] 
using the standard box size. In all the runs in Set D, the 
Reynolds and Coriolis numbers are kept above the critical 
values found by KMH11. In Runs D5, D4 and D3, we still 
observe the excitation of vortices, but in the rest, other types 
of large-scale flows are generated, suppressing the instabil- 
ity. Such large-scale flows are normally referred to as ba- 
nana cells, seen both in Cartesian (Chan 2001; Kapyla et 
al. 2004) and spherical geometries (e.g. Brown et al. 2008; 
Kapyla etal. 201 la, 201 lb). 

The structures extend through the whole convection 
zone and even penetrate to the overshoot region, see Fig. [4] 
where we plot a radial-vertical (xz) slice of the azimuthal 
velocity U y from Run D3 with #=30°. The convection and 
also the vortex tube, are strongly affected by rotation, which 
forces the structures to become inclined with the axis of ro- 
tation. 
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(lower panel) across the anticyclone for Runs Al and A2 
with differing box size. 




1 10 100 

fc/fc, 



Fig. 3 Kinetic energy spectra for the Runs Al (black 
lines) and A2 (red lines). Spectra taken at early times, when 
no vortices are yet excited, are plotted with dashed linestyle, 
and spectra from the vortex-state with solid linestyle. The 
wavenumber range of Run Al is scaled to match the one of 
Run A2. 
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Fig. 4 Two-dimensional slice, in the xz-plane, of the az- 
imuthal velocity U y from the run D3 at co-latitude 30 de- 
grees. The anticyclone generated in the run shows up as a 
structure spanning through the entire vertical extent of the 
computational domain, inclined with the rotation vector. 

The growth rate of the instability is reduced when the 
co-latitude is increased; therefore, it has been extremely dif- 
ficult to follow the Runs D3 and D4 up to saturation. The 
run at the pole (D5) having the largest growth rate has been 
run to a state very near saturation, when the growth of the 
horizontal velocity components slows down. In an earlier 
state, this run also exhibited both an anti-cyclonic and cy- 
clonic vortex; in the later stages, however, only the cyclone 
persists. The other runs still show both types of vortices, 
but this may still change as these runs are relatively further 
away from the saturated state in comparison to D5. 

3.3 Dependence on stratification 

We have made an attempt to quantify the effect of increas- 
ing stratification on the vortex instability by running one 
model where the stratification is almost 30 times larger than 
in our standard cases. Due to this, larger resolution in the ra- 
dial direction is needed; the amount of gridpoints has been 
increased from nz= 128 to 192, making these computations 
even more demanding than the rest included in this study. 
This run is identified as Run CI in Table [1] and can be ob- 
served to have a much higher Reynolds number (roughly 
94) than any other run, but clearly a lower Coriolis number 
(roughly 4.3) than the rest of the runs. Nevertheless, these 
numbers exceed the critical values found in KMH11, and 
therefore the vortex-instability is to be expected, unless the 
increased stratification has a significant effect on its excita- 
tion conditions. 

Indeed, we observe the instability in Run CI, gener- 
ating a stronger cyclonic and a weaker anti-cyclonic vor- 
tex, although the growth rate of it is reduced in compar- 
ison to the other runs due to the lower Coriolis number. 
We have been able to calculate this model up to roughly 
2000 turnover times, but it is clear that the system is still 
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Fig. 5 Azimuthal cuts through the computational domain 
at a; = —1.1 for various depths. The bluer the color, the 
nearer the bottom of the convection zone the cut is taken. 
Red colors represent cuts near z m , and the black line a cut 
taken exactly at z m . The yellow colors represent cuts taken 
near the top of the convection zone. 

far from saturation. Nevertheless, the temperature anomaly 
measured at this relatively early stage is already compara- 
ble (maximally close to 5%) to the standard box calcula- 
tions very near saturation. The vortices, again, occur at the 
very largest scales of the box, and their size does not vary 
significantly as a function of depth, even though the sys- 
tem is more strongly stratified. The temperature anomaly, 
on the contrary, varies monotonically through the convec- 
tion zone, see Fig. [5] where we plot a temperature cut in 
the y-direction through a cyclonic, cooler, region located at 
x = —1.1 for various depths. From this figure it is evident 
that the temperature contrast is monotonically increasing as 
function of depth. In the models with weaker stratification, 
such an effect is not clearly visible. 

4 Conclusions 

In an earlier study (KMH1 1) we reported on the excitation 
of large-scale vortices in Cartesian hydrodynamical convec- 
tion models subject to rapid enough rotation. In that study, 
the conditions of the onset of the instability were investi- 
gated in terms of the Reynolds and Coriolis numbers in 
models located at the stellar North pole. In this study, we 
extend our investigation to varying domain sizes, increasing 
stratification and place the box at different latitudes. 

The effect of the increasing box size is to increase the 
sizes of the generated structures, so that the principal vor- 
tex always fills roughly half of the computational domain. 
The instability becomes stronger in the sense that the tem- 
perature anomaly and change in the radial velocity are ob- 
served to be enhanced. Also the growth rate, measured from 
the time evolution of the horizontal velocity components, is 
more than doubled when the boxsize is doubled. The model 
with the smallest box size is found to be stable against the 



instability even though the critical Reynolds and Coriolis 
numbers found in the earlier study KMH11 are clearly ex- 
ceeded, suggesting that a sufficient scale separation between 
the convective eddies and the smallest wavenumber of the 
domain is required for the instability to work. 

The instability can be seen up to a co-latitude of 30 de- 
grees, but at higher co-latitudes, the flow becomes dom- 
inated by large-scale flows known as banana cells. Such 
flows have earlier been found from Cartesian and spheri- 
cal convection models. The vortices are seen to align with 
the rotation vector, forming structures tilted but coherent 
through the entire convection zone, extending even to the 
overshoot region. The growth rate of the instability is ob- 
served to become lower with increasing co-latitude. 

Only very little variance of the temperature contrast 
across the vortices can be seen as function of depth when 
stratification is small. The instability can also be seen in a 
model with larger stratification. Unlike the weakly strati- 
fied cases, the temperature anomaly caused by the vortex 
sturctures is seen to depend on depth. 
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